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ABSTRACT 


A procedure is presented for utilizing a bi-grid stability analysis as a practical tool for predicting 
multigrid performance in a range of numerical methods for solving Euler and Navier-Stokes equa- 
tions. Model problems based on the convection, diffusion and Burger’s equation are used to illustrate 
the superiority of the bi— grid analysis as a predictive tool for multigrid performance in comparison to 
the smoothing factor derived from conventional von Neumann analysis. For the Euler equations, 
bi-grid analysis is presented for three upwind difference based factorizations, namely Spatial, Ei- 
genvalue and Combination splits, and two central difference based factorizations, namely LU and 
ADI methods. In the former, both the Steger-Warming and van Leer flux-vector splitting methods 
are considered. For the Navier-Stokes equations, only the Beam-Warming (ADI) central difference 
scheme is considered. In each case, estimates of multigrid convergence rates from the bi-grid analy- 
sis are compared to smoothing factors obtained from single-grid stability analysis. Effects of grid 
aspect ratio and flow skewness are examined. Both predictions are compared with practical multi- 
grid convergence rates for 2-D Euler and Navier-Stokes solutions based on the Beam-Warming 


central scheme. 



1. Introduction 

Multiple grids were first proposed in the form of two-grid level schemes to accelerate the conver- 
gence of iterative procedures by researchers like Federenko [1]. Full multiple grid methods were 
later introduced by Federenko [2] to solve the Poisson equation and the approach was generalized by 
Bakhalov [3] to any second-order elliptic operator with continuous coefficients. According to Stu- 
ben and Trottenberg [4], Hackbush in [5] also independently developed some fundamental elements 
of the multigrid method. Perhaps the most influential work on the application of multigrid methods 
to elliptic type problems is that of Brandt [6] who also proposed the use of local mode analysis to 
determine the smoothing rate of multigrid schemes. 

In local mode analysis, the spectral radius of a particular relaxation technique computed over only 
the high-frequency modes is used as a measure of the relaxation’s effectiveness in a multigrid 
scheme since, in this case, the role of relaxation is not to reduce the total error but to smoothen it out 
i.e. remove the high-frequency components. It is assumed that the high-frequency modes have short 
wavelength that are spatially decoupled and that all high-frequency waves are completely ’killed’ 
on the fine grid and are not visible to the coarse grids. This, however, is not always the case since the 
inter— grid processes also influence the convergence rate. Brandt [7] presented theoretical consider- 
ations for including the transfer processes in the local mode analysis in what is called the bi-grid 
method. Also, some theoretical background is given by Stuben and Trottenberg [4] on how to com- 
pute a more realistic amplification factor for multigrid methods based on the bi-grid analysis, where 
some convergence norms were computed for the Poisson and Helmholz equations. 

A number of works exist where the smoothing factor has been used to predict multigrid performance 
in practice. However, the bi-grid analysis is becoming more attractive because of its better accuracy 
and reliability, van Asselt [8] used the bi-grid analysis to determine the proper amount of artificial 
viscosity to add at different level of coarse grids in a multigrid application. Mulder [9], [10] has also 
used the bi-grid method to construct an effective semi-coarsening in a multigrid method that can 
solve the problem of strong alignment which often occurs in convection problems. To select a relax- 
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ation scheme for a multigrid method suitable for a parallel solution of a time-dependent problem, 
Horton and Vandewall [11] employed this technique using the heat equation as their model problem. 
The cause of the poor multigrid convergence rate that is experienced in high-Reynolds number 
flows, where the coarse grid corrections fail to approximate the fine grid problem well enough for 
certain components, has also been investigated by Brandt and Yavneh [12] using the bi-grid method. 
In an effort to develop an effective multigrid algorithm for Navier-Stokes solutions on an unstruc- 
tured grid with 0(N) complexity, Morano [13], and Morano and Dervieux [14] have used the bi- 
grid analysis on a 1-D model scalar convection equation with periodic boundary conditions. More 
recently, Ibraheem and Demuren [15] also presented some convergence norms for the Burger’s 
equation based on bi-grid analysis. 

Implicit numerical schemes are becoming very popular, since they allow large time steps for advanc- 
ing the solution of Euler and Navier-Stokes equations to steady state. However, only few works exist 
to show the effectiveness of multigrid methods especially when approximate factorization is 
introduced. Yoon [ 1 6] and Caughey [17], for example, used the smoothing factor and scalar convec- 
tion equation as a model for the Euler equations to investigate multigrid performance. Anderson et. 
al. [18], and Demuren and Ibraheem [19] have also computed the smoothing factors on the actual 
coupled Euler equations for some popular approximate factorizations. The latter work investigated 
the Navier-Stokes equations as well. 

The objective of the present work is to present a procedure for utilizing the bi-grid amplification 
factor as a more practical tool for predicting multigrid performance in a range of numerical methods. 
Bi-grid analysis, based on the von-Neumann type method, is first presented for 1-D convection and 
diffusion model problems, and the linearized Burger’s equation. Numerical results from practical 
multigrid solution of these problems are compared to both predictions from bi-grid analysis and 
smoothing factors derived from the more usual single grid analysis. Both analyses and practical 
computations are based on the following different time-stepping methods: Euler forward explicit 
scheme, Runge-Kutta multistage scheme, a fully implicit scheme, and the semi-implicit scheme. 
The influence of the Peclet number on the convergence characteristics of the different schemes is 
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investigated using the Burger’s equation. Finally, for more practical situations, multigrid perfor- 
mance of various approximate factorizations for the 3-D Euler and Navier-Stokes equations are ex- 
amined using the bi-grid stability analysis. For the Euler equations, bi-grid analysis is presented for 
three upwind difference-based factorizations and several central difference-based factorizations. In 
the upwind factorizations, both the flux-vector splitting methods of Steger-Warming and van Leer 
are considered. The central schemes include the Lower and Upper (LU) and ADI factorizations. The 
time-stepping algorithm for the Navier-Stokes equations is based on the Beam-Warming central 
difference scheme only. Practical multigrid solutions from numerical experiment on the ADI meth- 
od are also compared to both predictions from bi— grid analysis and smoothing factors. 

2. Bi-grid Analysis 

Consider a given differential problem which can be written as: 

L{u(x)} = fix) ; for x in Q (1) 

where L is a linear operator. A typical 2-level multigrid cycle solution to this problem will involve 
the following steps: 

(1) pre-relaxation on a fine grid using any technique Si, v 1 times 

(2) computation of the defect R 

(3) restriction of the defect to the coarser grid 

(4) exact solution of the error equation on the coarse grid 

(5) prolongation of the error onto and the correction on the fine grid 

(6) post-relaxation on the fine grid using any technique S 2 , times 

These can be represented for any intermediate solution w, by using usual operators as follows: 
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( 2 ) 


(1) w n+ 5 = 

(2) R=f-L h w n+ t 

(3) I*R 

(4) v H = L^irfR) 

(5) l h H v H + 

(6) w n+l = S%(l h h f v„ + w n+ 5) 
Combining these steps, we can write: 


w n + 1 = S^U^h ll hV ~ L h‘ y i' vv ") + ( 3) 

The steady-state solution (u) is not changed by the coarse grid correction scheme, thus 

u n+ 1 = S v *[I h H L£ l l”(f - L h S\ l u n ) + S\'u"] (4) 

Subtracting (2) from (1) and noting that e n+1 = u n+l - w n+1 gives: 

= (5) 

= K S\'e n 
= Me n 

where K = I 

m - y/a - iiJ-s'’ ?W' (6) 


Mis the bi-grid amplification matrix and K is the coarse grid correction matrix. It can be shown [4] 
that when linear operators are used for the restriction, 1%, and the prolongation, /#, transfer pro- 
cesses, the coarse grid correction matrix is not a convergent iteration matrix, i.e., 

g(K) = q(L- I h H ^ l Ih L h) * 1 (7) 

Hence, the fine grid smoothing steps Si, and S 2 are important for a convergent scheme. The spectral 
radius of the bi-grid amplification matrix ( ) and its l 2 norm can be used to predict the per- 
formance of a multigrid method. While the spectral radius measures the asymptotic convergence rate 
of the multigrid method, the l 2 norm measures the actual error reduction per iteration, ^ 

defined as follows: 
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X max_bg ~ max{e[J#(0)]} (8) 

j#(0) is the Fourier representation of the matrix M. A brief comment about © will be in order. Due to 
aliasing process, low-frequency modes will couple with the coarse grid Fourier modes and, thus, for 
any 0 1 = {6 x ,9 y ,0 z } such that - n/2 < 6 X , 6 y , 6 Z < n/2 there exists a corresponding set of har- 
monics up to an integer multiple of 2ji. For 1— D, 2 — D and 3 — D problems, we define © as the follow- 
ing set : 


1- D © = ± *)} 

2- D © = {(0*,0 y ), [e x ,e y ± x),($ x ± a, o y ),(e x ± n,e y ± w)} (9) 

3- D © = {(d x ,Oy,e z ),(e x ,0y,6 z ± 7t),(0 x ,6y ± 7t,0 z ),(0 x ,dy ±7l,e z ± 7t), 

[0 X ± 7l,0y,0 Z ),(Ox ± 7t,0y,d z ± 7t), [0 X ± 71,9 y ± 71, 0 Z ), 

( 0 X ± 71,0 y ± 71,0 z ± 7r)j 

Or more generally, 

d-D © = [e l ,e 2 ,e 3 e ld ) (10) 

(where d is the dimensionality of the space, and 0 1 , <9 2 , ... . & ld are permuted in a similar manner 
with the ± signs chosen such that the harmonics lie in the high-frequency range). 

Hence, based on the 0 components and on the number of degrees of freedom of the problem, q, 

M(0) is a 2 d q x 2 d q matrix. Thus, it is a 2x2 matrix for a 1-D scalar problem and 40x40 matrix for 
the Euler or Navier-Stokes equations in 3-D. The Fourier representation for the corresponding op- 
erators viz: smoothing factor, fine grid problem, interpolation, restriction and the coarse grid prob- 
lem can be constructed as follows [7]: 

, £(<9 2 ),... . ^6> 2 ^)] 2 d q x 2 d q 

2 d q X 2 d q 
2 d q X q 
q X 2 d q 

qXq 


$ = (fcf, $*') = dl 
l h = diag[£(0 i y l(© 2 ), ^(® 2</ )] 

^ 0l )’ ^M 02 1 


fh = 

l H 


(ID 


% - ifM] 

K = A 201 ) 
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The difference operator, LrflO 1 ), on the coarse grid is only qxq since the coarse grid problem is 
solved exactly. 

$ and t depend on the choice of the smoother and the governing equations, respectively. The transfer 
processes, however, are less problem-dependent Following [7], the Fourier symbol of the pro- 
longation operator based on a / ,A -order polynomial is given by: 

d 

l h H {O m )ld = Y\Vl ( cos m ~ 1> 2< * ( 12 ) 

i = 1 

where \p 2 {g) = ( 14 - £)/ 2 , V> 4 (£) = (2 + 3£ - £ 2 )/4, etc. are the 2nd and 4 th order interpolation 
functions, and <5 w is the Kronecker delta. We restricted our analysis to the 2nd order since it is more 
commonly used. The restriction operator is expressed as: 

2<S/H(0 m ) = l / I h H (0 m )] T * (13) 

T* in the above equation represents the conjugate transpose. The restriction operator is often the 
adjoint of the prolongation operator in practice. In this study, the corresponding full weighting is 
used for the restriction operation for the Euler and Navier-Stokes equations, while simple injection 
is employed for the model problems. In the latter case, the Fourier symbol for the restriction operator 
is simply unity. 

A description of how the Fourier representation M(0) can be constructed is given later for certain 
problems. 

3.Model Equations 

The model equations used in the present study are the conservation equations for the convection of a 
scalar, the diffusion of a scalar, and the linearized Burger’s equation which is essentially a convec- 
tion-diffusion equation. Each of these equations is integrated in time using (i) Euler forward-explic- 
it scheme, (ii) a Runge-Kutta multistage scheme, (iii) a fully implicit scheme and (iv) a 
semi-implicit scheme. 
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The model equations for convection, diffusion, and the linear Burger’s equation can be written as : 


convection : u t + cu x = 0 

diffusion : u t - = 0 ( 14 ) 

(convection-diffusion) Burger’s : u* + u a u* = vu]^ 

In the Burger’s equation, u 0 = constant is assumed in our analysis. Thus, it can be put in the fol- 
lowing non-dimensional form: 

Wf + U x = ~p^ Uxx (15) 

where Pe in the above equation is the Peclet number defined as follows: 

pe = ( 16 ) 

(D is an appropriate length scale) 

(i) Euler forward-explicit scheme 

The Euler explicit method can be applied to the above equations to yield the following general dis- 
crete form: 

u? +l = u? - AtR n (17) 

where R n represents the residual expressed as follows: 
convection : R n = —{u 1 - - m"_j) 

diffusion : R n = - “ 2u l + u ?-i ) < 18 > 

Burger’s : R n = ^(k? - w?^) - ~ 2m ? + u 1-i> 

Space discretization in the above formulations is based on first-order upwind differences for con- 
vection, second— order central differences for diffusion, and the corresponding combination in the 
Burger’s equation. First-order upwind differencing of the convective flux introduces inaccuracy 
due to too much numerical diffusion which may be of the same order of the natural diffusion in the 
Burger’s equation. If second-order central differencing is used for the convective flux, a second-or- 
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der accurate scheme can be obtained but with severe limitations on the Peclet number due to disper- 
sion errors. Although the addition of artificial viscosity could dampen the high-frequency 
oscillations at high Peclet numbers, it is highly problem dependent. A better approach to achieve a 
second— order accuracy while sustaining a smooth solution at the vicinity of shock or high gradients 
is to discretize the convective flux using higher-order upwind schemes, preferably in conjunction 
with some limiter. Hence, with a third-order discretization of the convective flux, a second-order 
accurate scheme for the Burger’s equation can be obtained with R n given by: 


(19) 


R " = 2 ~ <-■> ” _ 3 “ ? + 3 “ ? - 1 ” “ ? - 2> 

\ (w?-ui — 2tt? + M? ,) 

Ax 2 Pe ,+1 ' ,_1 

(ii) Runge-Kutta Multistage scheme 

With each of the above schemes integrated in time using the Euler forward explicit method, the time 
step was limited to a small range by stability considerations, thus making it inefficient for steady- 
state computations. A Runge-Kutta (RK) method was introduced by Jameson et. al. [18] to permit 
larger time steps to be taken. For an m-stage scheme, the time integration can be written as follows: 


u °i = 

Ilf = M; - 
M «+l _ u m 


k = l,m 


( 20 ) 


Note that with m = 1, the RK scheme reduces to the Euler forward explicit scheme and hence is 
sometime called RK1. Coefficients a k are optimized such that larger time steps can be used for 

faster convergence. 

Three different sets of coefficients for a 4-stage Runge-Kutta scheme are investigated in this study, 
in line with the earlier work of Morano [16]. These are the standard coefficients (RK4-S, 
a 1 =. 25, a 2 = • 3333, a 3 =.5 ,a 4 = 1), and the optimized coefficients of Lallemand (RK4-L, 
a, =. 11, a 2 =. 2766, a 3 =. 5 ,a 4 = 1) and van Leer (RK4-VL, a x =. 0833, a 2 =. 2069, 
a 3 =. 4265, a 4 = 1 ) . 
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(Hi) Implicit scheme 

An implicit time integration scheme in delta form can easily be formulated for each of our model 
problems. For example, the corresponding implicit formulation for the Burger’s equation with first- 
order accuracy is written as follows: 


’ - + ?&)}■■?-. + [> + ' it + 


*" -*?+«*-,) 


f27j 


zlx 

Jm? = U? + 1 — M? 


in the above formulation is called the implicitness factor. /? = 1 . 0 gives a fully implicit scheme. 


(iv) Semi— implicit scheme 

Iff} = 0 . 5 in equation (21) above we have a semi-implicit scheme. This reduces to the Crank-Ni- 
colson scheme if the overall spatial differencing is second-order accurate. 

Fourier Symbols 

For illustration, the bi-grid amplification matrix M(0) is constructed for the convection problem 
using the Euler-forward explicit scheme for the smoother 

Consider the discrete form of the operator L and let the step-by-step solution be characterized by 
Fourier modes (with periodic boundary conditions) 


u n = Uj. n e 6Ji 


( 22 ) 


Then each of the operators that forms matrix M(@) becomes: 
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$(<9 m ) = |l - + ^f[cos(0 m ) - 7sin((9 m )] 

L^O™) = - -^[1 - cos(<9 m ) + I sin(0 m )] 

l h ^e m ) = ^[1 + cos(6> m )] m = 1,2 

I^{Q m ) = 1 for injection 

L w = “ ^[l " cos^# 1 ) + /sin^# 1 )] 

where 0 1 = 0 X and O 1 = 6 X + x 


Thus, from Eq. (8) M(0) can be written as: 


A _[st&) o] p: n ^ 2 1p(0') on 
M(&) ~ [o S( 0 2 )J [K 2l *J[o S(0 2 )\ 

Kn = i- &e')i*(e%(e')/L H 

K n = - fae')%(e%{e 2 )/L H 
K lt = - 4 ( 0 2 )^ , ^( 0 , )/4 
Jfa = 1 - 4<9 2 )/*(0 2 )4(0 2 )/4 


(24) 


Note that 1 H is evaluated only at the fundamental frequency , hence it is lxl. The result obtained 
above is similar to that derived by Morano [16], although our presentation is more general and is 
more easily extended to multi-dimensions. 

Multigrid Implementation 

A simple two-level multigrid (V cycle) method was implemented to test the relative accuracy of the 
bi— grid amplification factor and the smoothing factor in predicting multigrid performance. The 
two-level algorithm consists of the steps given in section 2 and is recursively expressible as follows: 


ll 



Proc Multigrid ,/?",£) 

^ ^either u n+l = L^ x R n 


or 


u 


n + 1 


= S°V 


else 


u 


n + l «_ 51 u n 


(25) 


R n <- I%{R n ~ Lu n ) 
Multigrid (0, u H ,R n , k — 1) 


f n + 1 


, 71+1 


endif} 


+ IfjUn 


In the above, L and S stand for the discrete operator and relaxation scheme corresponding to each of 
the model equations and numerical schemes discussed in previous sections. For this two-level V 
cycle multigrid implementation, the exact solution of the residual equation is employed. Only one 
pre-relaxation with no post-relaxation is performed on the fine grid. 


Local Relaxation 

Bi-grid analysis is exact for problems with periodic boundary conditions since it is based on the 
Fourier method. However, the asymptotic convergence rate for certain multigrid solutions deterio- 
rates from the bi-grid prediction due to singularities such as discontinuity in material and /or solu- 
tions, and also due to the type and coefficients of the boundary conditions. Poor multigrid 
performance results since such singularities lead to too large a correction from the coarse grids in the 
localized region. To improve the performance of a multigrid solution, further relaxation can be per- 
formed on fine grid in the region of the singularities after applying the coarse grid correction. This 
local Relaxation is, infact, an extra post-relaxation but confined to only certain nodal points and car- 
ried out a few number of times. The extra computational work is negligible if only a few partial 
sweeps is involved. The convection dominated problems subject to Dirichlet boundary conditions 
that are considered here undergo high changes in gradient in order to satisfy the exit boundary condi- 
tions. Therefore, multigrid performance in these problems deviates from the results predicted by bi- 
grid. However, a few passes on fine grid over the boundary conditions and over the interior equation 
in some small neighborhood of the boundary (about 3 nodal points at the exit) is found sufficient to 
improve multigrid performance to the exact value predicted by bi— grid analysis. 
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Numerical Experiments 

Bi-grid amplification factor ( bg ), the smoothing factor ( ) and the practical asymptotic 

convergence rate ( Q mg ) of the multigrid scheme were obtained for the following test problems: 


(1) Convection problem with periodic boundary conditions, viz.: 


«(0,f) = ii( 1, f) 


u(x, 0) = sin2jix 


(26) 


(2) Convection problem with Dirichlet boundary conditions, viz.: 


«(<U) = 1 


«(1, 0 = 0 for t>0 ; u(x, 0) = sin2jrx 


(27) 


(3) Diffusion problem with similar Dirichlet boundary conditions as in (2) above 

(4) Burger’s equation with similar Dirichlet boundary conditions as in (2) above. 

The bi-grid amplification factor is obtained from Eq. (8) and the smoothing factor is obtained from 
the usual single grid amplification factor over the high frequency range n/2 < (9 1 ^ Jr as 

Xp sg = maxjpl^© 1 )] ] . In each case, sixteen Fourier modes are selected, and the associated ei- 
genvalues are solved for using linear algebra routines such as found in the IMSL library. The asymp- 
totic convergence rate of the multigrid experiments, on the other hand, is computed from [19]: 


” III *■' II 


. n2-nl 


(28) 


where|| i? nl || and || R n2 || are the l 2 norm of the residuals at time levels nl and n2, respectively. 
The pseudotime At to advance the convection and the diffusion problems to steady state is com- 
puted from CFL = ~ and d = 2 , respectively. CFL is the Courant-Friedrichs-Lewy num- 


ber and d is the diffusion number. For the Burger’s equation, Jr is computed from: 
At — min(oJjc , oA^Pe) 


(29) 
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where o is an appropriate parameter chosen to reduce to the diffusion number d at low P e number 
and to the CFL number at high Pe number. This choice ensures that the appropriate time step is used 


in each flow regime. Ax is computed from D/20. Preliminary tests showed that the same results are 
obtained with 40 or 80 points. 

The exact steady-state solution for the Burger’s equation, subject to the boundary condition type 
discussed above, is given by: 


u = u(0,t) 


1 - exp [Pe (tj - l)] 
1 — exp(— Pe) 


(30) 


It is valid for all range of Pe considered in this study. 

Results for the Model Equations 

Figures 1 and 2 show results of the analyses of the 1-D convection equation using the Euler forward 
explicit scheme. The model problem of Fig. 1 has periodic boundary conditions whereas that of Fig. 
2 has Dirichlet boundary conditions. The bi-grid analysis gives perfect prediction of practical multi- 
grid performance in the former, whereas the smoothing factors from the single grid analysis are 
much too high. Both methods of analysis ignore boundary effects, so the same predictions are ob- 
tained in Figs. 1 and 2, and the analyses predictions are strictly correct only for problems with period- 
ic boundary conditions. This is confirmed in Fig. 2(b) where the asymptotic multigrid convergence 
rate is now much worse than predicted by the bi-grid analysis. The reason for the degradation of the 
multigrid performance is the singularity which appears near the exit in Fig. 2(a). This degradation in 
performance could be cured with a few local relaxation sweeps [15], as shown in Fig. 2(c). Each 
sweep had marginal computational cost and 5 sweeps were sufficient to bring the multigrid perfor- 
mance for the Dirichlet problem in line with that with periodic boundary conditions and the predic- 
tion of the bi-grid analysis. Clearly the Euler forward explicit scheme does not have good 
convergence properties except for CFL numbers close to 0.5, and it is divergent for CFL numbers 
greater than 1. Better convergence properties are achieved with Runge-Kutta (RK) schemes. Three 
4-stage RK schemes were analyzed, and the results are shown in Fig. 3 for the 1-D convection prob- 
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lem with periodic boundary conditions. With optimized coefficients Fig. 3(c), convergence could be 
obtained for CFL numbers up to 3. Further, bi-grid amplification factors below 0.4 are obtained for 
the range of CFL numbers from 0.5 to 2.5. There is also perfect agreement between the results of the 
bi-grid analysis and the practical multigrid convergence rates. Similar multigrid results were ob- 
tained by Morano [16]. Fig. 4 shows the result for the Dirichlet boundary conditions. In this case the 
multigrid convergence rates at higher CFL numbers are much better than predicted by either method. 
Clearly, the boundary effects are stronger with the RK scheme and there is no simple way to account 
for them in the analyses. Figure 5 shows results for a fully implicit scheme and for the semi— implicit 
Crank Nicolson scheme, for the 1-D convection equation. Although both schemes are stable for the 
whole range of CFL numbers, the Crank Nicolson scheme suffers from very poor convergence rate at 
high CFL numbers. 

Results for the 1-D diffusion equation are presented in Figs. 6-8. Dirichlet boundary conditions are 
applied throughout, and the steady state solution is shown in Fig. 6(a). In each case, the bi-grid anal- 
ysis gives perfect agreement with the multigrid convergence rate whereas the smoothing rate ob- 
tained from the single grid analysis is consistently too optimistic. On the whole, the predicted 
convergence rates for each method are similar to corresponding one obtained from the convection 
equation if the diffusion number, d is replaced by the CFL number in the latter. Clearly, if the goal is 
to achieve rapid convergence to the steady state, the fully implicit scheme with high d or CFL num- 
ber is the obvious choice. 

The linearized Burger’s equation represents a mixed convection-diffusion problem. The whole 
range of model type from pure diffusion to pure convection can be obtained simply by varying the 
Peclet number from a very small value to a very large value. Computed results for 4 values of Pe 
( 10 - 4 , 20, 100, 10 6 ) are presented in Figs. 9-12, for the various discretization schemes considered 
here. The exact solution at the steady state is shown in Fig. 9(a), for the Dirichlet boundary condi- 
tions u(0,t) = 1, u(l,t) = 0. For high values of Pe, there is a singularity near x=l. As explained in 
section 6 local relaxation is performed to reduce the adverse effect of this singularity on the overall 
multigrid convergence rate. The results for the first— and second— order Euler time explicit schemes 
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are presented in Figs. 9 and 10. In each case the bi-grid analysis gives quite good prediction of the 
multigrid convergence rate. On the other hand, single— grid analysis gives too optimistic estimates at 
low Pe and too pessimistic estimates at high Pe. The second-order scheme shows much poorer con- 
vergence rates, especially at high Pe. The results for the fully— implicit and semi— implicit schemes 
ate presented in Figs. 1 1 and 12. The superiority of the fully— implicit scheme is confirmed, especial- 
ly for high Pe flows. For <3 (or CFL number) greater than 10, it is close to a direct solver 

with A -»• 0. In these cases too, the bi-grid analysis agrees quite well with the practical multigrid 
convergence rate, except near a = 1 in the semi-implicit scheme at high Pe. Because of the lim- 
ited range of a where the convergence rate is much less that 1, the semi-implicit Crank Nicolson 
scheme is not a viable method for obtaining steady solutions for the model problem. If the main inter- 
est is rapid convergence to steady state, then the fully-implicit scheme at high values of o (or CFL 
number) will be optimum. 

Presendy bi-grid stability analysis has been presented for typical explicit and implicit solution 
methods for model problems which range from the diffusion equation to the convection equation and 
including the convection-diffusion equation at different Peclet numbers. For large scale practical 
computations, interest is really in solving the system of Euler or Navier— Stokes equations. In the 
following sections, the bi— grid stability analysis of fully-implicit schemes of Euler and Navier— 
Stokes equations are examined under various approximate factorization methods. 

4.Euler and Navier-Stokes Equations 

In order to extend the bi-grid analysis to the coupled equations of fluid flows, a discrete analog of 
these equations is formulated based on different approximate factorizations. The ADI factorization 
is formulated for the Navier-Stokes equations with the Euler equations as a degenerate case. Three 
different upwind factorizations and one central LU factorization formulated in [19] are, also, consid- 
ered. 

The 3-D Navier-Stokes equations in Cartesian coordinates can be written as 
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( 31 ) 


dQ 3(£ - Ey) 3(F - Fy) 3(G - Gy) _ q 
at + a* dy dz 


where is the solution vector and E, F, G are the conserved inviscid fluxes: 
Q = [g, gu, gv, gw, gef 

T 

E - [gu, gu 2 + p, guv, guw, (ge + p)u] 

F = [g>v, gvu, gv 2 + p, g>vw, (ge + p)v\ 

G = [g>w, gwu, pwv, gw 2 + p, (ge + p)w] 


(32) 


The viscous fluxes E Vy F v , G v are: 

E v = [o, |//(2 u x - v y - Wz), p(u y + v x ), p(u z + w*), T 

pv(u y + v x ) + pw(u z + w x ) + jfiu(2u x - v y - Wz) + fcTxj 

Fy = fa p(Uy + V X ), |//(2 Vy ~ U X ~ Wj, fX(V z + Wy), j. 

pu(u y + V*) + jHW(V Z + Wy) + |//v(2vy — u x - Wz) + fcTyj (33) 

Gy = |o, fi(W x + Hz), MVz + Wy), | J U(2W Z - Vy - Hj[), j. 

pu(w x + H z ) + //v(v z + Wy) + ^//w(2 Wj — Vy — H*) + £T z j 

In above, T = p/[gc v (y - 1)], and /> = (y - l)jpe - 0 . 5(u 2 + v 2 + w 2 )j. Also, Stokes hy- 
pothesis^ = - (2/3)//) has been assumed. With £ v , F v> G v set to zero, we recover the Euler equa- 
tions. 

Using the Beam-Warming scheme, the viscous fluxes are split directionally [20]. Following the ap- 
proach presented in Anderson et al. [21] for 2— D Navier— Stokes equations, analysis yields the fol- 
lowing ADI approximate factorization for the 3-D Navier-Stokes equations. Here, Euler time 
integration and constant fluid properties are assumed. 

[I + At0A - <5«K)][l + AKdyB - dyyS)](l + At(d z C - d^JAQ = (34) 

- At[Ad x - Rdxx ~ R^yx - R^ZX + Bdy - Sjd^y - 5<5yy - S-yd z y + C6 Z ~ ~ Fjdyj - Yd a ]Q 
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where the Jacobians A, B, C are dE/dQ, dF/dQ, BG/dQ, respectively. The analytical expression 
for the viscous fluxes are given in Deniuren and Ibraheeni [19]. The right-hand side resulted from 
linearization and from assuming the flux Jacobians to be locally constant. To damp the high-fre- 
quency waves that will arise due to central differencing, second— order implicit 
(D* = - s^tAxdxx) and fourth-order explicit (D e x = - t ( AtAx , b XX3 ^ artificial dissipations are 
added as diagonal matrix coefficients in the numerical examples. Thus, with similar dissipations 
added in the y and z directions Eq. (34) becomes 

[i + A t(d x A — dxxR — £^dxdxr)][l + At(dyB — dyyS — £^dydyy)]^I + A t(d z C — d a Y — SjAzd^^AQ 
= — A$Ad x — Ed xx ~ E{dy X — R 2 dzx + Bd y — S^dxy ~ Sdyy — S 2 dzy ( 35 ) 

+ Cd z — — Y 2 dyz ~~ Ydiz + e e (Ax^d xxxx + Ay^dyyyy + Az^d zzzz )jQ 

The corresponding factorization for the Euler equations becomes apparent if the viscous flux Jaco- 
bians/?, R v R 2 , S, S v S 2 , Y, Fj, Y 2 are set to zero. 

Other approximate factorizations that are considered in this work are those formulated for Euler 
equations in [18] and [19], viz: 

[I + At(6xA + + <5+A-)][l + At(dyB + + d, + 2T )Jl + At(KC + + <5 Z + C")]A<2 = ~ AtR " ( 36 > 


[I + At{6x A + + 6yB + + <5 2 ~C + )][l + At(d?A~ + <5 +5" + <5 Z + C-)]A(2 - - AtR n (37) 

[I + At0xA + + d +A- + dfC + )][ I + At(6yB + + d+ZT + d+C~j\AQ = - Atl R n (38) 

where R n = d x E + + 6}E~ + d y 'F + + d?F~ + d z G + + d + G~ (39) 

Jl + At0 x Ai + 6y Bi + d z Cj) + x^Atidx + d y + d z )] 

X [i + At(d x + A 2 + 6 y + B 2 + d z + C 2 ) - XiAtid? + dy + <5 Z + )] (40) 


SS — At[d x E + dyF + <5 Z G) “ XjAtiAX^dxXXX AjPdyyyy + A z d ZZZZ 'j 
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Eqs. (36), (37) and (38) arc upwind schemes, and are referred to as spatial, eigenvalue and combina- 
tion factorizations, respectively. The flux-vector splitting methods of Steger-W arming [22] and van 
Leer [23] are assumed. Eq. (40) is the Lower and Upper (LU) factorization. Here, the fluxes devised 
by Jameson and Turkel [24], viz: Aj = (A 4- IAI)/2 and A 2 = (A — IAI)/2, are used to achieve 

diagonal dominance. d + and d~ denotes forward and backward difference operators, respective- 
ly, and x 2 and x 4 are the artificial dissipation coefficients. 

Fourier Symbols 

The bi-grid amplification matrix M(0) is constructed from M = I — ItfLjj * ^ or 

ease of presentation, the Euler equations alone are selected for illustration, with the ADI central 
scheme used as the smoother. In this case, viscous fluxes/?, /? 2 » & ^ 2 , Y v F 2 areset 

to zero. The components operators of matrix M(@) are: 

A 

(i) The fine/coarse grid Operator L 
The Euler equivalent form of Eq. (31) is: 

§=-(f + f + f) + dissipmion (41) 


Thus, in quasi-linear form: 

«*» - - (''f + ”w * c $) + + 


(42) 


Holding A, B, C locally constant and employing second-order central differencing, the Fourier sym- 
bol of the fine grid problem on equal mesh size in all directions becomes: 


L A (<9 m ) = - -f[A sin(0f) + B sin((9^) + Csin(<9£>] + ^[cos(6>p + cos(<9£) + cos(0f) - 3] 


Ax 


16 Se 

Ax 


sin4 (^) + sin4 (x) + sin4 (x) m=l,8 


(43) 
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Note that 0 ” represent the f^ h element of the © m component (see Eq. (9-11)). 

For any arbitrary mode, Eq. (43) is a 40 X 40 matrix since each Jacobian is a 5x5 matrix and there are 
8 harmonics including the fundamental mode. The coarse grid problem is assumed to be a version of 
the original problem on the fine grid and the coarse grid is formed simply by deleting every other fme 
grid point. Thus, the mesh size and Fourier modes are { 2Ax, 20 1 } and its Fourier signature can be 
written as: 

Lf/2& 1) = - ^ [A sin(20 x ) + 5sin(20 y ) + Csin(20*)] + ^[cos(20,) + cos(2 6 y ) 

+ 005(20*)- 3] -^(sm 4 0* + sin 4 0 y + sin 4 0 z ) (44) 

In the above equation, only the fundamental mode, 0 1 = {Q Xj Oy, 0 Z } > is employed since the coarse 
grid problem is assumed to be solved exactly. Hence, this is only a 5 X 5 matrix. 

(ii) The relaxation Operator $ 

Each of the equations (35), (36), (37), (38) and (40) can be expressed as 

NAQ n - - L - — AtR n (45) 

von Neumann stability analysis is used on this system of linear equations by letting the step-by-step 
solution be characterized by 

Q n = (46) 

where A is the single grid amplification factor. Thus, Eq. (45) reduces to a complex generalized ei- 
genvalue problem of the form 

A A AAA / *<n\ 

Kx = XNx where K = N — L (47) 

The Fourier symbols of Nand'L, for our particular example, can easily be shown to be 
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1 JZJtffi / 

IX0" 1 ) = ^l{A sin(0f) + B sin(6>p + Csin(6>£)) + ^ I sin 4 -^- + sin 4 -^- + sin 4 -^- 1 (49) 


The Fourier symbols corresponding to the other approximate factorizations are documented in De- 
muren and Ibraheem [25]. For each harmonic, O m (m = 1, 8), Eq. (47) is solved to give five eigen- 
values from which the elements of $(0) are constructed. For example, if the eigenvalues 
corresponding to the mode @ l = {0j(,0y,0z} are A = {A l ,A 2 ,A 3 ,A 4 ,A 5 }, then, from Eq. (11), 
§(0t) = yiL The effective fine grid smoothing operation is obtained by raising the smoothing ma- 
trices to the power of v 1 andv 2 , the pre- and post-smoothing counts, respectively. 


(iii) The Transfer Operators /^and/f 

For a second-order interpolation, the Fourier symbol of the prolongation operator, from Eq. ( 1 2), is: 
7^(<9 m ) = |[l + cos(0j*)][l + cos(<9£)][l + cos^p] (50) 


The restriction operator, ??,is computed from this equation andEq. (13) assuming full-weighing. 

Based on the above operators, M(0) is assembled from M = S^(L — I .A symbolic 

form is given in Appendix A. It is an 8x8 block matrix of which each elemental block is a 5x5 matrix. 

Solution Procedure 

The eigenvalues for the bi-grid matrix i$(0) are computed from Eq. (8) over fixed Fourier modes 
to obtain the amplification factor. Sixteen modes are selected, in the range — n/T. <, © l < nfl. 
The smoothing factor is also computed from the generalized eigenvalue problem (47) over only the 
high-frequency modes ;r/4 < I0 1 ! < Ji/2d&Xp_ sg = max(IAI). In each case, the eigenvalues are 
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solved for using the linear algebra routines such as found in the IMSL library. Uniform flow is as- 
sumed with Mao =0.8, zero yaw (a y ) and angle of attack ( a a ), andy =1.4. Further, the grid 
spacing is assumed to be uniform in all directions. Effects of aspect ratio and flow skewness are also 
investigated. The time-step and Reynolds number are calculated from 


At = 


CFL 


jii ■ m , 


M 


Tz 


-I- C 



Jy 2 Jz 2 


(51) 


+ Ay ® + 

Some other pertinent definitions used are as follows: 

I VI = 7m 2 + v 2 + w 2 , Moo = Jjp , v = u tan(a>,) , w = u tan(a a ) 




(52) 


(53) 


Practical multigrid solutions are obtained for Euler and Viscous flows around a circular cylinder us- 
ing the PROTEUS computer code developed at NASA Lewis Research Center. The FAS-FMG (full 
approximate storage-full Multigrid) algorithm applicable to non-linear systems of equations was 
implemented in the PROTEUS code [26]. Based on two levels, the asymptotic convergence rate of 
these flows were computed from Eq. (28). The Reynolds number based on the cylinder diameter is 
20 and the Mach number is 0.2. The grid size for the Euler flow is 25x49 and for the Viscous flow is 
49x49. In each case, the grid was clustered such that the aspect ratio varied from 1.5 to 3.8 for the 
Euler flow and 0.5 to 12.2 for the Viscous flow. Further, the pre- and post-smoothing counts are 1 
and 0, respectively, as are assumed in the analyses too. 

Results for the Euler and Navier-Stokes Equations 

Figure 13 shows the convergence results for the 3-D Euler equations using the upwind schemes. The 
computed values for the smoothing factor (A^_$j) and bi-grid amplification factor (A max ^)for the 
spatial, eigenvalue and combination factorizations based on the Steger-Warming flux-vector split- 
ting are shown in figures 13(a), 13(b) and 13(c), respectively. Both factors predict instability for the 
spatial split, especially for CFL number beyond 5.0. In the eigenvalue and combination factoriza- 
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tions, better convergence characteristics are observed although the smoothing factor’s prediction is 
slightly more optimistic. For these two factorizations, bi-grid analysis predicts near instability at 
CFL number above 25 whereas the smoothing factor predicts unconditional stability for all CFL 
numbers. Figures 13(d)— 13(f) show predictions for multigrid performance of each factorization us- 
ing the van Leer flux-vector splitting. Except for the spatial factorization, all the schemes are pre- 
dicted unconditionally stable for all CFL numbers by both bi-grid and smoothing factors. The 
spatial factorization is stable only for CFL numbers below 1 2 and possesses better convergence char- 
acteristics at CFL number below 8 than the other two factorizations. From both analyses, i.e. from 
(Xftjg) and (A max bg \ van Leer flux-vector splitting gives better convergence characteristics than 
the Steger— Warming method for multigrid procedures. It is observed that the present results of the 
smoothing factors for the van Leer method are similar to those presented by Anderson et. al. [ 1 8] , and 
Demuren and Ibraheem [19]. 

Results for the 3-D Euler equations using the LU approximate factorization with central difference 
approximation and various levels of second- and fourth-order artificial viscosities, x 2 and x 4 , are 
shown in Figs. 14(a)— 14(c). Without the addition of second— order dissipations, i.e. x 2 — 0, the co- 
efficients x 4 = 0 . 3 yields the optimal results (see fig. 14(a)). From figures 14(b) and 14(c), bi- 
grid and smoothing factors predict that an appropriate combination of x 2 and x 4 (especially when 
x 4 > x 2 ) can significandy improve the performance of the LU scheme when used as a relaxation 
scheme for multigrid. Also for all levels of dissipation, the smoothing factors estimates are more 
optimistic than the bi-grid results especially at lower CFL numbers. 

The convergence characteristics for the 3-D Euler and Navier-Stokes equations for different levels 
of artificial dissipation and Reynolds numbers are shown in figures 14(d)— 14(0 and 15, using the 
Beam-Warming (ADI) central difference scheme as the baseline solution algorithm. With no dis- 
sipation added to the Euler equations (fig. 14(d)), the bi-grid analysis predicts instability for all CFL 
numbers, while the smoothing factor predicts stability for CFL numbers below 15. From figures 
14(e) and 14(f), optimal multigrid performance is predicted by the bi-grid analysis for dissipation 


23 



levels of t e = 0.5 and e- = 1.0. These results are similar to those obtained for the Navier- 
Stokes equations at Re=10® (see figs. 15(d)— 15(f)). With Reynolds number of 100 and no dissipa- 
tion, both bi-grid and smoothing factors predict stability for certain range of CFL numbers although 
the latter is more optimistic. Also at this Reynolds number, the optimal dissipation levels are 
s e = 0 . 5 and e f = 1.0. 

All computations have been based on zero yaw and angle of attack, and also on uniform grid spacing 
in all directions. Sensitivities of convergence characteristics to flow skewness and aspect ratio are 
studied using the ADI central scheme at Reynolds number of 100, and dissipation levels of 
= 0 . 5 ande, = 1.0 .The results are shown in figures 16 and 17. Generally, convergence char- 
acteristics are improved with increases in yaw angle at zero angle of attack although the range of 
stable CFL numbers becomes smaller (figs. 16(a)- 16(c)). From figures 16(d)— 16(f), no significant 
difference is observed in the convergence results when the yaw and angle of attack are set equal to 
each other. From figure 17, the convergence characteristics become worse with increases in grid as- 
pect ratio. 

In order to ascertain the suitability of bi-grid and smoothing factors in predicting multigrid perfor- 
mance in complex flows, asymptotic convergence rate are computed from practical multigrid solu- 
tions of 2-D Euler and viscous flows around a circular cylinder. The steady-state solutions for these 
flows are shown in fig. 18(a). Rather than evaluating the corresponding bi-grid and smoothing fac- 
tors from uniform flow conditions, as in previous cases, they are computed at each point in the flow 
field, thereby accounting for the variation in flow properties. Figs. 18(b) and 18(c) show estimates 
from both analyses based on the computed frozen coefficients of the Euler and viscous flows, respec- 
tively. These results are also summarized in Table I, and are compared with the asymptotic conver- 
gence rate measured from the practical multigrid computations. For both flow problems, smoothing 
factor deviates more from the practical solution than does the bi-grid factor. 
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Concluding Remarks 

Bi-grid stability analysis has been presented for typical explicit and implicit solution methods for 
model problems which range from the diffusion equation to the convection equation and including 
the convection-diffusion equation at different Peclet numbers. Bi-grid amplification factors were 
compared with smoothing factors and multigrid convergence rates. The predicted bi-grid amplifica- 
tion factors agree quite well with the asymptotic convergence rate of the multigrid method. The 
smoothing rate of the relaxation scheme obtained from a local mode analysis on a single grid is not an 
accurate predictor of the multigrid convergence rate. For multigrid performance in large scale practi- 
cal computations, bi-grid amplification factor and smoothing factor were computed from the system 
of 3-D Euler and Navier-Stokes equations. Various approximate factorization methods that are pop- 
ular in practice are considered. In typical practical multigrid solutions of 2-D Euler and viscous flow 
problems, bi-grid analysis is found to give a better prediction of the convergence rate than the 
smoothing factor obtained from a single grid analysis. 
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TABLE I: Convergence characteristics of 2-D Euler and Viscous flows around a cylinder 



Euler 

Viscous flow 

CFL 



Qmg 



Qmg 

0.5 

0.88 j 

0.94 

0.99 

0.95 

0.96 

0.99 

1.0 

0.80 

0.92 

0.98 

0.91 

0.94 

0.98 

2.0 

0.76 

0.91 

0.96 

0.85 

0.93 

0.96 

4.0 

0.81 

0.90 

0.93 

0.77 

0.92 

0.94 

6.0 

0.84 

0.90 

0.92 

0.76 

0.91 

0.93 

8.0 

0.87 

0.90 

0.92 

0.81 

0.91 

0.92 

10.0 

0.89 

0.90 

0.94 

0.84 

0.91 

0.92 

12.0 

0.91 

0.90 

0.95 

0.87 

0.91 

0.92 
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APPENDIX A 


The Bi-grid Amplification Matrix A?(0) 




The diagonal elements are: 


Mn = I 

M* = I - 1 

m 33 = i - Wki^^i0^ 2 i0^ 


- I - /^0 8 f/r(<9 8 ^(0 8 )^ , (6» 8 ^^ 1 
and the off-diagonal elements are: 

where, for example, 


m* - - iy,eii(()'Ue'$(e=)%{et;' 

m„ = - 

m u -- 

m„ - - faeisnet^i^st; 1 

Each element is a 5x5 matrix corresponding to the 5 dependent variables in equation (31). 
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Fig. 2: 1-D Convection Equation (a) Steady solution (b) Convergence Characteristics 
without local relaxation (c) Convergence Characteristics with local relaxation 
(Euler forward explicit; DIrichlet B.C’s) 
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(a) RK4— S 




Fig. 3 : Convergence Characteristics for 1— D Convection E qua t i on , 
(a) Standard (b) LeDamad (c) van Leer coefficients 
(4— Stage Runge Kutta; Periodic B.C*s) 
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Fig. 4: Convergence Characteristics for 1— D Convection Equation 
(4-Stage Runge Kutta; Dirichlet B.C*s; van Leer coefficients) 
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(b) 



Fig. 5 : Convergence Characteristics for 1 — D Convection Equa t i on 
(a) Implicit (b) Semi-implicit time integrations 
(Wrichlet B.C*«) 
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Fig. 7: Convergence Characteristics for 1— D Diffusion Equation 

(4-Stage Runge Kutta; Dirichlet B.C*s; van Lear coefficients) 
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d 

Fig. 8 : Convergence Characteristics far 1— D Diffusion Equation 
(a) Implicit (b) Semi-implicit time integrations 
(Dirichlet B.C*s) 


39 


















0.0 0.2 0.4 0.6 0.8 1.0 


0.0 0.2 0.4 0.6 0.8 1.0 



Fig. 10 : 1— D Linear Burger’s Equation (a)-(d) Convergence Characteristics 
(Euler forward explicit; 2nd 0 accurate) 
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(a) Pe = 10 


(b) Pe - 20 



(c) Pe = 100 


(d) Pe - 10® 



Fig. 11 : i-D linear Burger** Equation (a)-(d) Convergence Characteristics 

(ImpHdt time integration) 
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(a) Pe = 10' 


(b) Pe = 20 




Q 


(c) Pe - 100 



(d) Pe = 10 # 



a 


Fig. 12 : 1— D Linear Burger’s Equation (a)-(d) Convergence Characteristics 
(Semi— implicit time integration) 
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(a) Spatial Factorization 
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0.4 
0.2 
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Ctl 

(b) Eigenvalue Factorization 


0.0 

0 10 _ 20 30 

cfl 

(d) Spatial Factorization 

1.2 
1.0 
0.8 
X 0.6 
0.4 
0.2 
0.0 

0 10 - 20 30 

cfl 

(e) Eigenvalue Factorization 





Steger & Warming 


1.2 
1.0 
0.8 
X 0.6 
0.4 
0.2 
0.0 

0 10_ 20 30 

cfl 

(f) Combination Factorization 
van Leer 



jpjg ^3 3_i) Euler Equation* uiin| upwind schemes (u)~(f) Convergence chamctaristic* 

(v*-! ; v*— 0) 
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(c) <2=3, <4=2 (f) e,=l. e<=2 

LU Factorization ADI Factorization 


pig. 14 3-D Euler Equations using central schemes (a)— (f) Convergence characteristics 

(^■1 ; v* ""O) 
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(c) e.=l, e i =2 (f) «,-l. £ i= 2 

Re = 10 2 Re “ 10 * 


pjg. 15 3_d Naviar-Sloka* Equation* using central schemes (a)— (f) Convergence characteristics 

(v>-l ; ^-0) 
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(c) 0^=60°, <x a =0° (f) 0^=60°. a a =60° 
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jjg is 3 — d Navier— : Stokes Equations using central schemes (a)— (I) Convergence characteristics 

(Re-100, S.-0.5 , Cj-1.0, v 1 -! , 0) 


47 

















'Ay=Ax/Az=10 



(c) Ax/Ay= 1000, Ax/Az= 1 (f) Ax/Ay=Ax/Az= 1000 

Aspect Ratio Aspect Ratio 


Pig 17 3-D Navier— Stokes Equations using central schemes (a)-(Q Convergence characteristics 

(Re-100, c,-0.5, C|-L0. i^-t v*-0) 

48 














0 5 10 15 20 0 5 10 15 20 

cfl cfl 

(b) Inviscid flow (c) Viscous flow 


Fig. 18 2— D Euler end Navier-Stokes flows around a circular cylinder using ADI scheme 

(^-iy-0. C.“l,c t =2) 
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